all_p2b <- readRDS("temp/precinct_block.rds")

pop <- function(s){
  print(s)
  library(tidycensus)
  library(data.table)
  library(dplyr)
  if(!file.exists(paste0("temp/block_pop_2010_", s, ".rds"))){
    
    cs <- get_decennial(geography = "county", variable = "P001001", state = s,
                       year = 2010, sumfile = "sf1")

    s_pop <- rbindlist(lapply(substring(cs$GEOID, 3, 5), function(c){
      get_decennial(geography = "block", variable = "P001001", state = s,
                    year = 2010, sumfile = "sf1", county = c)
    }))
    saveRDS(s_pop, paste0("temp/block_pop_2010_", s, ".rds"))
  }else{
    s_pop <- readRDS(paste0("temp/block_pop_2010_", s, ".rds"))
  }
  return(s_pop)
}


cl <- makeCluster(8)  
registerDoParallel(cl)

clusterExport(cl, list("pop"))

slist <- unique(filter(fips_codes, state_code <= "56")$state)

j <- rbindlist(parLapply(cl, slist, fun = pop))

all_p2b <- left_join(all_p2b, select(j, GEOID, pop = value), by = c("GEOID10" = "GEOID"))

bgs_level <- all_p2b |> 
  mutate(GEOID = substring(GEOID10, 1, 12)) |> 
  group_by(state, precinct, GEOID, year) |> 
  summarize(pop = sum(pop, na.rm = T))

c16 <- readRDS("temp/census_bgs_16.rds") %>%
  mutate(year = 2016) %>%
  select(GEOID, year, latino, asian, nh_white, nh_black, median_income, median_age,
         pop_dens, some_college)


c20 <- readRDS("temp/census_bgs_19.rds") %>%
  mutate(year = 2020) %>%
  select(GEOID, year, latino, asian, nh_white, nh_black, median_income, median_age,
         pop_dens, some_college)

census <- bind_rows(c16, c20)

pct_level <- left_join(bgs_level, census)

pct_level <- pct_level |> 
  group_by(state, precinct, year) |> 
  summarize(across(c(latino, asian, nh_white, nh_black,
                     median_income, median_age, pop_dens, some_college),
                   ~ weighted.mean(., pop, na.rm = T)))

saveRDS(filter(pct_level, !is.na(precinct)), "temp/pct_demos.rds")
#################################################
# read coded killings keep only well-coded killings
full_set <- readRDS("temp/geocoded_shootings_from_gva.rds") |> 
  mutate(date = as.Date(date),
         year = floor(year(date) / 2) * 2) |> 
  filter(home + dv == 0,
         year %% 4 == 0) |> 
  select(id2 = id, year, date)

precincts <- bind_rows(
  readRDS("temp/shootings_precincts_2016.rds"),
  readRDS("temp/shootings_precincts_2020.rds")) |> 
  mutate(GEOID = toupper(GEOID)) |> 
  rename(precinct = GEOID)

full_set <- left_join(full_set, precincts)

full_set <- left_join(full_set, readRDS("temp/pct_demos.rds")) |> 
  mutate(pre = (year == "2016" & date <= "2016-11-08") | (year == "2020" & date <= "2020-11-03"))

full_set <- as.data.table(full_set)

set_pre <- full_set[full_set$pre,
                    .SD[date %in% max(date)], by=list(year, precinct)]

set_pre <- set_pre[, .SD[1], by = .(year, precinct)] %>% 
  mutate(treated = T)

set_post <- full_set[!full_set$pre,
                     .SD[date %in% min(date)], by=list(year, precinct)]

set_post <- set_post[!(paste0(set_post$precinct, set_post$year) %in% paste0(set_pre$precinct, set_pre$year)),
                     .SD[1], by = .(year, precinct)] %>% 
  mutate(treated = F)

full_treat <- bind_rows(set_pre, set_post) %>% 
  select(precinct, id = id2, date, year,
         median_income, nh_white, nh_black, median_age, pop_dens, latino, asian,
         some_college, treated, share_dem) %>% 
  mutate(d2 = ifelse(year == "2020", as.integer(date - as.Date("2020-11-03")),
                     as.integer(date - as.Date("2016-11-08"))),
         t16 = year == "2016")

l <- rdrobust(y = full_treat$share_dem, x = full_treat$d2, p = 1, c = 0, cluster = full_treat$id,
              covs = select(full_treat,
                            latino, nh_white, asian,
                            nh_black, median_income, median_age,
                            pop_dens, some_college, t16)
              )
## create naive RDIT plot for Figure 4
j <- rdplot(y = full_treat$share_dem, x = full_treat$d2, c = 0, p = 1,
            covs = select(full_treat,
                          latino, nh_white, asian,
                          nh_black, median_income, median_age,
                          pop_dens, some_college, t16)
            )[["rdplot"]]
binned_dots_data <- j[["layers"]][[1]][["data"]] |> 
  select(x = rdplot_mean_bin, y = rdplot_mean_y)

left_line <- j[["layers"]][[2]][["data"]] |> 
  select(x = x_plot_l, y = y_hat_l)

right_line <- j[["layers"]][[3]][["data"]] |> 
  select(x = x_plot_r, y = y_hat_r)

all_dots = select(full_treat, x = d2, y = share_dem)

j[["labels"]] <- j[["labels"]][-1]
j <- j +
  geom_point(aes(x = d2, y = share_dem), data = full_treat, shape = 21,
             color = "transparent", fill = "gray", alpha = 0.5) +
  theme_bc(base_family = "LM Roman 10") +
  labs(x = "Days between mass shooting and election",
       y = "Democratic vote share") +
  scale_y_continuous(labels = percent) +
  scale_x_continuous(breaks = c(-180, -90, 0, 90, 180),
                     labels = c("Shooting occurs\n180 days\nbefore election", "90",
                                "Shooting occurs on\nelection day",
                                "90", "Shooting occurs\n180 days\nafter election"),
                     limits = c(-190, 190))
j

h <- j[["layers"]][[5]]

j[["layers"]][[5]] <- j[["layers"]][[1]]

j[["layers"]][[1]] <- h
j
saveRDS(j, "temp/rd_plot_demshare.rds")

for(i in outpaths){
  ggsave(paste0(i, "rd_plot_demshare.png"), height = 4, width = 6, units = "in")
}

fwrite(binned_dots_data, "data_for_figures/figure_4_blue_dots.csv")
fwrite(all_dots, "data_for_figures/figure_4_gray_dots.csv")
fwrite(left_line, "data_for_figures/figure_4_line_left_of_zero.csv")
fwrite(right_line, "data_for_figures/figure_4_line_right_of_zero.csv")

